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Abstract 

Grid staggering is a well known remedy for the problem of velocity/ 
pressure coupling in incompressible flow calculations. Numerous incon- 
veniences occur, however, when staggered grids are implemented, partic- 
ularly when a general-purpose code, capable of handling irregular three- 
dimensional domains, is sought. In several non-staggered grid numerical 
procedures proposed in the literature, the velocity /pressure coupling is 
achieved by either pressure or velocity (momentum) averaging. This ap- 
proach is not convenient for simultaneous (block) solvers that are preferred 
when using multigrid methods. A new method is introduced in this pa- 
per that is based upon non-staggered grid formulation with a set of virtual 
cell face velocities used for pressure/velocity coupling. Instead of pressure 
or velocity averaging, a momentum balance at the cell face is used as a 
link between the momentum and mass balance constraints. The numeri- 
cal stencil is limited to 9 nodes (in 2D) or 27 nodes (in 3D) both during 
the smoothing and inter-grid transfer, which is a convenient feature when 
a block point solver is applied. The results for a lid-driven cavity and a 
cube in a lid-driven cavity are presented and compared to staggered grid 
calculations using the same multigrid algorithm. The method is shown to 
be stable and produce a smooth (wiggle-free) pressure field. 
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1 Introduction 

Multigrid methods are used in a number of applications in fluid dynamics, 
usually by applying the Full Approximation Scheme [1], Incompressible flow 
calculations usually employ a staggered grid because of its strong coupling 
between the pressure and the velocity field (e.g. [2j. For complex geometries, 
however, as well as For calculations in non-orthogonaT coordinates, the use 
of a staggered grid is a serious obstacle to efficient and well structured 
computer coding [3]. Additional complexities arise when a block-solver is 
used; for example, variables cannot be easily grouped into cell-bound blocks 
due to different node count. Some authors resort to asymmetric nodal 
clusters [5] while others update a symmetric block of variables around the 
cell centre node thereby updating face velocities twice in each relaxation 
sweep [5, 6]. Various levels of decoupled relaxation axe also common. These 
include distributive relaxation, where all momentum equations are solved 
together and the pressure field is solved separately [1, 7], and sequential 
schemes that update variables throughout the flow field one by one [8, 9]. 
Some comparative studies of block versus sequential relaxation give no clear 
preference [10, 11]. There is a greater consensus that grid staggering is a 
necessary burden, particularly in the context of multigrid methods ([12, 13, 
14, 15, 16] and even [1]). Comparison studies of staggered and non-staggered 
methods are sometimes conflicting in their assessment of the accuracy and 
stability of any given method. While some authors demonstrate that non- 
staggered methods match the staggered ones using both criteria ([13, 16, 
17]), others question it ([18]). Despite this, the majority of finite volume 
incompressible calculations use staggered grids. The main reason may be 
that existing non-staggered grids increase rather than lessen the complexity 
of the staggered grid calculations. For example, the method of Rhie and 
Chow [19] (adopted by [13, 14, 15]) requires that both the nodal and 
c e ll fa ce velocities are stored. Moreover, in a multigrid context, both the 
nodal and the face velocities need to be restricted ]8] , requiring even more 
computational work. Also, the computational cluster extends beyond either 
9 or 27 point stencil in two- or three-dimensional formulations respectively 
for the first order discretisation and even more if the higher order methods 
are used. 

The considerations mentioned above motivated the present contribution 
for a method suitable for block solvers on an irregular three-dimensional 
domain using a non-staggered grid. 


In this paper a brief description of the multigrid procedure based on 
a new non-staggered method is given in section 2 and the multigrid 
implementation in section 3; the test cases and results are presented in 
section 4, followed by the discussion section where relative merits of the 
method are assessed. 


2 The new non-staggered method 


A transport equation for a general set of transported variables u in the 
volume fi bounded by a boundary S can be expressed in an integral form 
suitable for finite volume formulation 
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where p is the density, u t is the velocity component in the direction and 
ra j is the component of unit normal to the boundary S. When (1) is applied 
to the momentum balance of a viscous incompressible fluid, the set u is a 
velocity vector u = {uj,j = 1, ..., d} (d being the problem dimension), with 
the corresponding diffusion coefficient T = p and the source terms 



( 2 ) 


in the absence of external volume forces. The extension to other trans- 
portable properties (such as enthalpy, mass fraction, etc.) is straightforward 
by augmenting the vector u to include new variables and defining appropri- 
ate source terms and the diffusion coefficients. In the following presentation 
a three-dimensional implementation will be used. 

The momentum equations are discretised using the hybrid (central/up- 
wind) difference scheme [20] although higher order schemes can be 
employed 1 . The pressure field is resolved by means of mass conservation for 
the control volume around the node in a symmetric block manner as used 
by Vanka [5] for the staggered grid, although the extension to the line block 
around the node in a symmetric block manner as used by Vanka [5] for 
the staggered grid, although the extension to the fine block formulation is 

J For a multigrid implementation of a second order upwind scheme on a staggered grid 
see e.g. [21]. 


straightforward. The estimation of the face velocities that are substituted 

into the mass conservation equation is obtaine d by dis cretizing the mo- - 

mentum equation over a half-length volume around each cell face, directly j 

involving the nodal pressure and velocities, while the lateral velocities are 
obtained by averaging values from the nearby nodes (see Fig. 1). More 
details on the coefficient generation are given in [22]. 

The principle of discretizing the face velocity using a half-size cell is 
applied also by Schneider and Raw [3], although in their method the 

coefficients of the face velocity are treated implicitly by incorporating them j 

into the nodal velocity coefficient matrix. To ensure positive definiteness of 
the nodal velocities coefficient matrix, Schneider and Raw had to truncate 
the momentum equation applied to the face velocities [3, 16]. In the 
present method, the face velocities are explicitly expressed in terms of 
the surrounding nodal values and used in the continuity equation for 
the pressure correction calculation. The implication of this step is that 
the family of face velocities satisfies both the momentum and the mass 
conservation exactly at the positions where the convection velocities in 
a general transport equation are required. On the other hand, as an 
average of the (tentative) nodal velocities they are readily available without 
requirement for a permanent storage. 

The boundari es of the flow field are coincident with the cell faces, 
enabling the definition of a set of boundary nodes there. This practice 
ensures consistency between the global mass balance of the w hole calculation 
domain and the local mass balance of each cell, but calls for special 
treatment if Neumann boundary conditions are to be used. If the zero- 
gradient condition for the normal velocity is discretised in a usual way 

d(xiUi) ( XiTii)p - (g tn^ fon* _ 

where subscripts b and inn denote boundary and the first inner node, 
respectively, the flow rate through the boundary will be linked to the 
velocity that does not belong to the mass preserving field, resulting in 
poor overall mass conservation. The correct way to implement Neumann 
boundary condition in this case is to use the face velocity. This way, the 
local and global mass balance become fully compatible. There is no need 
for any special treatment of the Dirichlet boundary conditions where the 
face velocities coincide with the boundary and are assumed known. 
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3 The multi-grid implementation 

In the multigrid context, the nonlinear equation (1) can be expressed as 

C( u) = F (4) 

by grouping all terms that will result in a coefficient matrix (within a 
Newton iteration cycle) into the operator C and the remainder into the 
source term F as in [13, 23]. A more common practice of expressing Eqn. (4) 
as homogeneous (by absorbing F into £(u)) [1, 24] is found by the present 
authors to be somewhat confusing, especially when defining residual transfer 
to the coarse grid. 

The discretised (sparse, positive definite) Eqn. 4 for the grid l is linearised 
by a Newton iteration [24] 

L l u l = F l (5) 

and relaxed by a block Gauss-Seidel method. 

The updates of the variable set 

u' = a(diag(L)) _1 R* (6) 

are expressed in terms of the residual R/ = F* — it u*, the inverse of the 
coefficient matrix diagonal (diag(A)) and the underrelaxation coefficient 
a. Variables at the node i,j, k are then updated by = u i,j,k + u 'i,j,k- 
Restriction is accomplished by grouping a cluster of eight cells into one. 
This leads to the following operator 

<f>I,J,K = + &+lj\Jfc +<f>i,j+l,k + </>i,j,fc+l + 

<f>i+lj+l,k + 0iJ+l,Jb+l + 0t+l,j,fc+l + 0i+ l^'+l.fc+l)) ( 7 ) 

where I = 2i — 2, . . . The same operator is applied both to variable and 
residual restriction. After both the variables and residuals are transferred 
to the next coarser grid (l — 1), Eqn. (5) is approximated as 

L Z-1 (u' _1 ) =F* - \ (8) 


F *- 1 


= F*- 1 - (Fjj- 1 - L q -1 ^ -1 ) + 7^ _1 R* 


where 
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is the equivalent source term on the coarse grid. The restriction at 
Neumann boundaries is carried out using a divided form of the boundary 
conditions [l] . For the velocity component perpendicular to the boundary, 
an additional correction is made to preserve the mass flow rate through the 
boundary. 

Prolongation is carried out by tri-linear interpolation using a seven point 
stencil, shown here for one cell and with injection only: 


<l>ij,k = g (3<f>I,J,K + + 4>I,J-1,K + 1) (!0) 

with i = (J + 2)/2, . . . The injection upon the first visit to the fine grid 
(FMG cycling is assumed) and the fine grid correction are done as 

“first = or u new = Uol d + P/— l(U{— l-U 0 ,Z-i). (11) 

4 Test cases 


The non-staggered method presented in this paper is compared with 



Gauss-Seidel algorithm of Vanka [5] . For both methods the coding and data 
structures are of the same style. 


The flow in a three-dimensional cavity with a moving top is used as a 
first test case. The residual norm history is shown in Fig. 2. The rate of 
convergence obtained when calculating on a staggered grid is comparable 
with the results of Vanka [10] where 12 work units (w.u.) were necessary for 
a two orders of magnitude residual reduction. In our calculations 14 w.u. 
was necessary for the staggered grid calculation and 18 w.u. for the non- 
staggered calculations. 1 However, the change in slope of the non-staggered 
residual may indicate that the full potential of multigrid acceleration has 
yet to be achieved. 

In a second test case, a cube is inserted in a cavity (Fig. 3), forcing the 
flow to negotiate this asymmetric three-dimensional obstacle, partly by the 
velocity magnitude change, partly by flow separation. It is believed that 
this flow geometry serves as a good test of the pressure/ velocity coupling 

1 Or 23 w.u. for the same residual decrease; however, this is more arbitrary, because 
of the much lower initial residual at the finest grid. 
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because the major force behind the flow adjustment is the pressure field. 
The residual history, (Fig. 4) indicates very similar convergence rates for 
the staggered and non-staggered calculations. The resulting flow field in a 
symmetry plane (Fig. 5) indicates well resolved separation bubbles around 
the cube comers. 


5 Discussion and conclusions 

The new method of incompressible flow calculation using non-staggered 
grids and its multigrid implementation are examined for suitability in a 
complex flow field geometry. The presence of two sets of velocity values, 
both of which satisfy the (discrete form of the) governing equations increases 
the overall level of accuracy for a given grid size, although this remains to 
be quantified. 

In the numerical experiments performed so far the method proved to be 
stable, without any tendency to produce an oscillating pressure field, which 
is a common feature of some non-staggered methods [18]. The method 
permits discretization on a trivially coarse grid (with a single node in the 
interior), which is very convenient in a FAS multigrid implementation, 
because it allows the coarse grid to contain the lowest number of nodes. 
Thus significantly coarser grids can be used in complex geometries. For 
example, in the case of a cube in a cavity (see the previous section) the 
coarsest grid (6x6x6 nodes) has only one control volume located between 
the cube and the cavity wall at one side. If the calculation method required 
two nodes at minimum, the overall node count at the coarse grid would 
increase eight times, thereby substantially increasing the work needed to 
obtain exact solution at the coarsest grid. 

Various tests performed so far always produced smooth solutions both 
in velocity and pressure, which indicate a high ellipticity measure of the 
proposed method. The analytical evaluation of the ellipticity measure 
remains to be carried out (following e.g. [25, 16]). 

The amount of computational work of the proposed method is slightly 
larger that of the Rhie and Chow [19] method. It is comparable to 
the work in the SCGS method of Vanka, requiring the same amount of 
work to calculate face velocities and pressure coefficients and, in addition, 
the calculation of the nodal velocity coefficients, i.e. approximately 25% 
increase in two-dimensional and 14% in three-dimensional calculations. This 



overhead exists only for the simplest flow problem because any additional 
variable that is solved permits nodal velocity coefficients to be reused (with 
proper scaling of the diffusion part). . 
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Figure 1: The layout of a non-staggered grid. Only the nodal variables 
require storage. 


3D Cavity Flow 
5 grids: 2x2x2 ... 32x32x32 





Figure 3 : The grid for a flow around the cube in a lid-driven cavity. 


Flow around the cube in a cavity 
4 grids: 6x6x6 ... 48x48x48 



Work unita 

Staggered grid Non -Staggered grid 

Figure Mass residual history for the flow around the cube in a lid-driven 
cavity. Cavity Re=400. 
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